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ABSTRACT 

How did Pluto's recently discovered minor moons form? Ward and Canup propose an elegant 
solution in which Nix and Hydra formed in the collision that produced Charon, then were caught into 
corotation resonances with Charon, and finally were transported to their current location as Charon 
migrated outwards. We show with numerical integrations that, if Charon's eccentricity is judiciously 
chosen, this scenario works beautifully for either Nix or Hydra. However, it cannot work for both 
Nix and Hydra simultaneously. To transport Nix, Charon's eccentricity must satisfy ec < 0.024; 
otherwise, the second order Lindblad resonance at 4:1 overlaps with the corotation resonance, leading 
to chaos. To transport Hydra, ec > 0.7i?pi uto /acharon > 0.04; otherwise migration would be faster 
than libration, and Hydra would slip out of resonance. These two restrictions conflict. Having ruled 
out this scenario, we suggest an alternative: that many small bodies were captured from the nebular 
disk, and they were responsible for forming, migrating and damping Nix and Hydra. If this is true, 
small moons could be common around large Kuiper belt objects. 
Subject headings: 



1. INTRODUCTION 

The recent discovery of Pluto's two minor moons Nix 
and Hydra (Weaver et al. 2006; Buie et al. 2006) presents 
interesting puzzles. We summarize current observational 
data in Table 1 . 

• The orbits of the moons are nearly circular and 
nearly coplanar with Charon's orbit. 

• Nix, the inner minor moon, lies just inward of the 
4:1 resonance with Charon, while Hydra is close to 
the 6:1 resonance. 

The Pluto-Charon system is doubly synchronized and 
circularized - it must have gone through significant tidal 
evolution in the past. In the currently favoured the- 
ory for the formation of Charon (Canup 2005; McK- 
innon 1989), a giant impact chipped off a piece of the 
proto-Pluto, leaving Charon on an eccentric orbit close 
to a rapidly spinning Pluto. Subsequent tidal evolution 
slowed down Pluto's spin, pushed out Charon to its cur- 
rent position, and damped its orbital eccentricity. This 
is similar to how Earth's Moon is thought to have formed 
and evolved. Tides on Pluto pushed Charon out to its 
current position in ~ 2 x 10 7 (Qp/100) yrs, where Qp 
is Pluto's tidal quality factor. Charon's eccentricity 
eventually decayed to zero on a comparable timescale, 
assuming that Charon's tidal parameter Qc is not too 
different from Qp. 4 By contrast, Nix and Hydra likely 
cannot evolve tidally at their current positions in the age 
of the Solar System (Stern et al. 2006; Lithwick & Wu 
2008). 
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3 Estimates of various tidal timescales are listed in Appendix B. 

4 Its eccentricity would initially have grown if Qc/Qp exceeds 
a number that is of order unity; otherwise, its eccentricity would 
have decreased monotonically. 



2. FORCED RESONANT MIGRATION (FRM) 
2.1. The FRM Scenario 

Ward & Canup (2006) propose an elegant scenario to 
account for Nix and Hydra's observed orbital proper- 
ties. In their scenario, the minor moons were formed 
as byproducts of the collision that formed Charon. Nix 
was then caught into Charon's 4:1 corotation resonance, 
and Hydra into the 6:1 corotation resonance. As Pluto's 
tides pushed out Charon, Nix and Hydra remained in 
resonance and so they too were forced to migrate to- 
wards their current orbits. In this scenario, Nix and 
Hydra must have been caught into the corotation res- 
onance, and not into any of the other sub-resonances at 
4:1 and 6:1, because migration in other sub-resonances 
would have excited the eccentricities of Nix and Hydra 
to values much larger than are observed. 

A corotation resonance can transport a particle only 
if Charon's eccentricity ec is sufficiently large, because 
the resonant libration time must be shorter than the 
time for the resonance to migrate a distance of order 
its width, and the libration time increases with decreas- 
ing ec whereas the width decreases. The more strin- 
gent constraint is set by Hydra, which requires ec > 
0.7(100/Q P ) 1 / 5 (i?p/a c ), where ac is Charon's semima- 
jor axis (Ward & Canup 2006). Therefore when Charon 
reached its current orbit at ac — 17-Rp, it must have still 
had an eccentricity of ec > 0.04. As ec was subsequently 
damped by tides, the width of the corotation resonances 
shrunk to zero, and Nix and Hydra escaped from res- 
onance. Such a history for Charon's orbit is plausible, 
given the uncertainties in the tidal parameters. 

Ward & Canup (2006) briefly address the question of 
how the minor moons were initially trapped into corota- 
tion resonances. If Nix and Hydra were produced in a 
collision, then their free eccentricities were initially large, 
and it would have been unlikely that they had just the 
right orbits to end in corotation. But if a lot of debris was 
produced in the collision, and if this debris was highly 
collisional, then it is possible that the debris settled into 
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TABLE 1 
Observed parameters 



[units] a 


Pluto 


Charon b 


Nix c 


Hydra c 


orbital period 


[To] 


1 


1 


3.89 


5.98 


semi-major axis 


[Rp] 


1.96 


16.81 


41.82 


55.65 


mass 


[M P ] 


1 


0.1165 


7.8 x 10~ 6 


1.8 x 10~ 5 


eccentricity 




0.0 


0.0 


0.002 


0.005 


inclination 


[deg] 


96.14 


96.14 


96.18 


96.36 



a Pluto's radius and mass are Rp = 1164 km (Young & Binzcl 
1994), and M P = 1.3 X 10 25 g. Orbital period of Pluto-Charon 
binary is T c = 6.3872 days. 

b Charon's radius and mass from Sicardy et al. (2006). 

c Nix and Hydra's radii and masses asssume a Charon-like albedo 
of 0.35 (Weaver et al. 2006) and a density of 2g/ cm 3 (Charon-like); 
masses will be ~ 20 times higher (and radii ~ 2.7 times larger) if 
albedo is 0.04 (comet- like). Orbital parameters are from Buie et al. 
(2006). 

a cold disk with very little free eccentricity, 5 and that 
Nix and Hydra's free eccentricities were damped by this 
disk. It is also possible that Nix and Hydra formed from 
this debris. 

We note that it is not possible that Nix and Hydra 
were formed outside of corotation, and then were colli- 
sionlcssly caught into corotation by a migrating Charon. 
The reason for this is that the contribution of a corota- 
tion resonance to the Hamiltonian is a cosine term with 
constant coefficient. 6 By the symmetry of such a term, 
the energy lost when the particle approaches the sepa- 
ratrix (the "balance of energy") must equal the energy 
gained when the particle leaves the separatrix, and cap- 
ture is impossible (e.g., equation 68 in Henrard 1982). 

A potential difficulty with the FRM scenario not ad- 
dressed by Ward & Canup (2006) is the effect of the 
forced secular eccentricity. If Nix and Hydra resided 
in corotation resonances when Charon's eccentricity was 
high (ec > 0.03), then their forced secular eccentric- 
ity was comparable to e C - Yet their current eccentric- 
ity is <C 0.03. How was this eccentricity damped? The 
answer is that their secular forced eccentricity tracks 
e c , and that as ec decays, so does the forced secular 
eccentricity — as long as the decay time of ec is longer 
than the secular precession time. To verify this, we start 
from the secular part of the Hamiltonian for a test par- 
ticle in the presence of Charon (see Appendix A) 



Z {Z) = -nn ( e c B 1 Z + Z + B 2 \Z\ 2 



(1) 



where Z = ee~ m is the particle's complex eccentricity, 
n is its orbital frequency, /x = Mc/Mp and Bi,B 2 are 
sums of Laplace coefficients. The evolution of Z is given 
by Hamilton's equation (eq. [A16]), 



dZ _ dH scc 
~dt~ 1 dZ* 

= —2in(j, I e c 



B,Z 



(2) 
(3) 



5 The particles try to settle into the "coldest" orbits possible, 
with nested orbits that do not intersect each other if possible. 

6 See, for example, the Co term in equation (A26). To be more 
precise, the coefficient of a corotation resonance does not depend on 
the eccentricity of the particle, i.e., Nix or Hydra. It does depend 
on the eccentricity of Charon, and on the semimajor axis through 
the Laplace coefficient. However, the coefficient changes very little 
as the resonance sweeps over the particle. 



The forced secular eccentricity is eforced = ~ecBi/2B 2 , 
and the secular precession frequency is lo scc = 2n[iB 2 . 
If ec decays exponentially on timescale r e , i.e., ec(t) — 
and if Z\ t= o — e IO rccd, then the solution of 



ec\t=oe 



-t/r e 



equation (3) is 



-t/T. 



Z = eforced 1 1=0 : — j ( 4 ) 

Therefore when t 3> r e and ec has decayed to 0, the 
particle's eccentricity will be negligibly small as long as 
^socTe 3> 1, proving our assertion. 

2.2. Numerical Simulation of FRM 

Figure 1 demonstrates numerically that the FRM sce- 
nario works for Hydra when Charon's migration trajec- 
tory is chosen judiciously. We start Charon on an orbit 
with ec = 0.1, ac = lQRp and force it to migrate out- 
ward to ap = 17 Rp with a = a/10 6 yrs, keeping the 
eccentricity constant. Then ec is forced to decay to zero 
in 10 6 yrs. A massless test particle representing Hydra is 
initially placed in the center of Charon's 6:1 corotation 
resonance. (See §3.2 for a description of how we start 
off a particle in the corotation resonance). The orbital 
motions of all bodies arc numerically integrated with the 
SWIFT package (Levison & Duncan 1994), using the hi- 
erarchical Jacobi symplectic integrator of Beust (2003). 
We further modify it to allow semi-major axis and ec- 
centricity evolution due to external forces, following the 
approach of Lee & Pealc (2002). 

Figure 1 shows that Charon indeed pushes out Hy- 
dra. Hydra remains within the corotation island until 
Charon's eccentricity falls to a very small value. Hydra's 
eccentricity tracks that of Charon's along the way (eq. 
[4]). Fig. 1 also shows that FRM depends extremely 
sensitively on the initial conditions of the test particle. 
A second particle started off at exactly the same location 
as the one above but with a velocity larger by 0.1% (and 
therefore falling outside the narrow corotation island) is 
quickly ejected during the migration of Charon. In fact, 
the width of the corotation island is so narrow that it 
is comparable to the size of Hydra at the start of our 
experiment. 

3. RULING OUT FRM 

FRM can migrate cither Hydra or Nix to their cur- 
rent orbits. But it cannot simultaneously migrate both 
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Fig. 1. — Numerical simulation illustrating the success of the 
Ward & Canup (2006) forced resonant migration scenario. Charon 
is forced to follow a simple migration trajectory: it is first pushed 
outward with a timescale of 10 6 years from 10 to 17 pluto radii 
(the left panel shows its semi-major axis as a function of time as a 
solid line), keeping its eccentricity fixed at 0.1 (solid line in right 
panel). After this we force its eccentricity to decay in 10 6 yrs. A 
test particle initially placed in Charon's 6:1 corotation resonance 
is resonantly migrated outward (dotted curves and open circles). 
Notice that its eccentricity falls off following that of Charon. The 
value of its secular forced eccentricity is 0.366c at the 6:1 location. 
A particle starting at the same location but having a velocity that is 
larger by 0.1% is quickly ejected (star symbols). To reduce scatter, 
we have plotted the test particle values only when all bodies have 
true anomaly A ~ 0. 

of them. To transport Hydra, Charon's eccentricity must 
be greater than a critical value; otherwise, Hydra would 
slip out of resonance. To transport Nix, Charon's eccen- 
tricity must be less than a second critical value; other- 
wise, the 4:1 corotation resonance would be destroyed by 
resonance overlap. These two constraints conflict. We 
discuss them in turn. 

3.1. Lower Limit on ec from Hydra's migration 

To be able to transport Hydra in resonance, the migra- 
tion time across the width of the 6:1 corotation resonance 
must be longer than the libration period in the resonance, 

6 2 ^c < ^ , (5) 

J lib 

where Aa^b is the resonance width, Xi;b is the libration 
period, and dc is Charon's tidal migration rate. 

We take the expressions for the libration width and li- 
bration period from equations (8.58) and (8.47) of Mur- 
ray & Dermott (1999), and use the Kaula formula (§6.3 of 
Murray & Dermott 1999) to obtain a resonance strength 
for the 6:1 corotation resonance of fd ~ 0.02 (as defined 
in eq. [8.32] of Murray & Dermott 1999). Adopting the 
tidal timescale for orbital expansion as in equation B6, 
the above transport condition translates to a lower limit 
for ec, 

Note the weak dependence on k 2 p (Pluto's tidal Love 
number) and Qp (Pluto's tidal quality factor). The 
above limit has also been given in Ward & Canup (2006). 
We have performed numerical experiments to confirm the 
numerical coefficient. We plot this limit in Figure 2. 

3.2. Upper Limit on ec from Resonance Overlap At 
Nix's Orbit 
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Fig. 2. — This figure shows why FRM cannot work. The solid 
lines show the lower limit of ec below which tidal migration of 
Charon is too rapid for forced resonant migration for Qp = 100. 
These lines are given by eq. [6] for the 6:1 resonance, and an anal- 
ogous expression for the 4:1 resonance. The dotted lines show the 
same limits, but for Qp = 1000. The dot-dashed horizontal lines 
show the maximum values of ec, above which individual corota- 
tion resonances are destroyed. For the 4:1 resonance one must have 
ec < 0.024. FRM is feasible for each resonance only if ec lies be- 
tween these two limits (shaded region, for Qp = 100). The lower 
limit for the 6:1 resonance lies above the upper limit for 4:1 reso- 
nance. Hence Charon could not resonantly migrate Nix and Hydra 
simultaneously. 

In this subsection, we demonstrate with two numeri- 
cal experiments that Charon's 4:1 corotation resonance 
exists only if ec < 0.024. Otherwise, resonance over- 
lap destroys the corotation resonance (as shown in §4). 
Figure 2 shows that the constraint ec < 0.024, together 
with the constraint from §3.1, rules out the FRM sce- 
nario: it is impossible to satisfy both of these constraints 
simultaneously. 

To demonstrate numerically the destruction of the 4:1 
corotation resonance, one first needs a robust method 
to find it when it exists. This is not entirely trivial be- 
cause the mass ratio Mc/Mp is not terribly small, and 
neither are Charon's and Nix's eccentricities ec and e. 
Hence an expansion of the disturbing function will not be 
very accurate. Instead, our method is based on finding 
the "coldest" orbits of test particles around the Pluto- 
Charon binary. A disk of infinitesimal particles that col- 
lide inelastically will naturally tend to the coldest or- 
bits, i.e. the orbits that minimize the velocity dispersion 
within the disk. What are these coldest orbits? Consider 
first the case that Pluto and Charon orbit each other on 
circular orbits. Then the coldest orbits are the periodic 
orbits, which are orbits that, in the rotating reference 
frame of the binary, close on themselves after a single 
loop around the binary. The generalization of periodic 
orbits to the case of non-zero ec are invariant loops (Ma- 
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ciejewski & Sparke 1997; Pichardo et al. 2005), which 
may be understood as follows. If we take a snapshot 
of the position of an orbiting particle every time Pluto- 
Charon reach some predefined orbital phase (periapse, 
say), then in general the snapshots would fill out a two- 
dimensional region in the r, 9 plane, where r, 9 are the 
particle's radius and azimuth. But the coldest orbits — 
the invariant loops — are those in which the snapshots 
trace out a one-dimensional closed loop. Far from strong 
resonances, neighbouring invariant loops do not intersect 
each other, or themselves. But near strong resonances, 
they often do intersect. 

To find the corotation resonance, we find the set of in- 
variant loops around an eccentric Pluto-Charon binary, 
largely following the procedure of Pichardo et al. (2005). 
We express the orbit of the binary in terms of the ec- 
centric anomaly, expanded to fourth order in ec- The 
motion of the test particle is evolved in the barycen- 
tric frame with a 4th-order Rungc-Kutta integrator with 
adaptive step-sizes (Press et al. 1992), with an error tol- 
erance of 10 -8 . A test particle is initially launched far 
from any strong resonances on the periastron axis (x- 
axis) when Charon is at periapse. We assume that its 
orbit is symmetric with respect to the cc-axis, so the only 
unknown is the velocity v y . This is initially guessed us- 
ing the local Keplerian value. The positions of the test 
particle at each subsequent periapse passage of Charon 
are recorded, for a total of 10 4 binary periods. These are 
separated into 2, 000 angular bins and the radials disper- 
sion within each angular bin is co-added. 7 We use the 
bisection technique (typically within a range 6% of the 
initial guess) to find the correct v y that minimizes this 
dispersion. The resulting 1-D curve r{9) is an invariant 
loop. We then proceed to find the next loop which is 
closer to the resonance, using the previously obtained v y 
as the initial guess. The closer the loops, the better the 
guess, and the more reliable the convergence. 

Fig. 3 depicts the invariant loops we obtain for a test 
particle near the 4:1 location. When ec < 0.024, 4:1 
corotation islands are clearly visible; while when ec = 
0.024 and beyond, all islands disappear. 

Fig. 4 shows the results from a second experiment 
that demonstrates the destruction of the 4:1 corotation 
resonance for large ec- In this experiment, we insert a 
test particle into the center of the 4:1 corotation reso- 
nance when ec is small, and then slowly raise the value 
of ec- As ec rises above 0.024, the motion of a parti- 
cle initially trapped inside the 4:1 corotation resonance 
becomes chaotic and the resonance angle circulates. 

Similar experiments performed for the 6:1 corotation 
find the limiting ec to be 0.24 in that case. 

4. RESONANCE OVERLAP DESTROYS THE 4:1 
COROTATION RESONANCE 

We seek a better understanding of how and why the 
4:1 corotation resonance is destroyed when ec > 0.024. 
We model the evolution of a test particle near the 4:1 
resonance of an eccentric Charon-Pluto binary with a 
truncated disturbing function, keeping only secular and 
4:1 resonant terms. The particle's equations of motion 
are compactly encoded by the Hamiltonian derived in 

7 This is different from Pichardo et al. (2005) in that they only 
use the dispersion near the x-axis. We find that our treatment 
gives a faster and more reliable convergence. 
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Fig. 3. — Invariant loops near the 4:1 location, showing the 4:1 
corotation resonance (as 4 distinct islands) when ec = 0.01 (left) 
and its absence when e c = 0.024 (right). These loops represent or- 
bits that are the "coldest" possible (no free eccentricity). Here, we 
record particle positions (barycentric radius in unit of Rp, position 
angle 9 as measured from Charon's periapsis) whenever Charon 
is at periapsis. The resulting surface-of-section is either a con- 
tinuous line (no resonance) or disjoint islands (in resonance). For 
ec = 0.01, other 4 : 1 subresonances lie beneath the corotation res- 
onance (not shown here); when ec = 0.024, subresonances overlap 
and there are no periodic orbits within the resonant region. 

Appendix A (eq. [A26]), which has two degrees of free- 
dom. The momentum and co-ordinate of the first de- 
gree of freedom {p a , q a } are related to the test particle's 
semimajor axis and mean longitude via equations (A13) 
and (A14). And the momentum and co-ordinate of the 
second degree of freedom are combined into the complex 
canonical variable z, which is the test particle's free com- 
plex eccentricity (i.e., the complex eccentricity after sub- 
tracting the forced secular eccentricity, eqs. [A15] and 
[A 19]). The equations of motion for these two degrees 
of freedom are just Hamilton's equations (eq. [A27]), 
which we shall numerically integrate. Since this is a time- 
independent Hamiltonian with two degrees of freedom, 
phase space may be mapped out with Poincare surfaces 
of section (e.g., Henon & Heiles 1964) 

We first consider the corotation term (the cq term in 
Hamiltonian A26) in isolation, discarding the Lindblad 
terms (the C1-C3 terms), which leaves 



H 



^c e c cosq a - fxB 2 \z\ 



(7) 



with the two degrees of freedom decoupled from each 
other; c and B 2 are order-unity constants whose values 
are listed in Appendix A. Hamilton's equations for q a ,p a 
yield q a — — 48p a and p a = —ficoe c smq a , the same as 
for a pendulum. Hence the angle q a will either librate 
around the center of the resonance (where q a = p a = 0), 
or, for large enough \p a \, it will circulate (Fig. 5, top 
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Fig. 4. — Maximum ec above which the 4:1 corotation reso- 
nance disappears. We place Charon at its current position, and a 
test particle at 4:1 corotation resonance. As ec is increased grad- 
ually (solid curve in top panel, e/e = 1/10 5 years), the barycentric 
eccentricity of the test particle rises even faster, until ec = 0.024 
(marked by the dotted line) after which resonance overlap destroys 
the corotation resonance — the resonant angle 4>a-.i = 4A — Ac — 3<I>c 
suddenly starts circulating. The particle is released from the reso- 
nance and displays chaotic motion. To reduce scatter, we only show 
test particle data when both bodies have true anomalies A ~ 0. The 
same break-down for the 6:1 resonance occurs at ec = 0.24. 

left panel). Hamilton's equation for z is z = —2i/iB2Z, 
showing that the free eccentricity has constant amplitude 
and a circulating phase (Fig 5, top right panel). 

As in §3.2, we seek the "coldest" orbits, i.e. the in- 
variant loops. An invariant loop appears as a fixed point 
in the z-plane surface of section. To see this, recall that 
an invariant loop is a ID curve of the test particle's r(9) 
(radius vs. azimuth) whenever Charon reaches periapse. 
When Charon reaches peripase, the test particle has 8 

9(t)=*& (8) 

r ^ = l + 8p a (i)-Re(e i0 W(^) + e forced )) , (9) 

valid to first order in z and p a , where a ros and e forcod 
are constants defined in the Appendix. Since the z-plane 
surface of section has 9 = q a = 0, we see that an invariant 
loop must be a fixed point in that section; otherwise, 
there would be an infinite number of different values of z 
(and hence of r) at 6 = 0, and the orbit would not trace 
out a ID curve in the r-6 plane. In summary, to find 
the invariant loops, we first find the fixed point in the z- 
plane surface of section. We then plot the r vs. 9 for this 

8 More generally, we should write 8 = q a /& + /2, where j is an 
integer. The four corotation islands in the r-9 plane are produced 
by different values of j . 



Fig. 5. — Orbits arising from the corotation resonance only, from 
numerical integrations of Hamiltonian (7)'s equations of motion 
with ec = 0.01. Top two panels show surfaces of section at a fixed 
value of the energy. Cross-hatched regions arc incompatible with 
the chosen energy and are forbidden. Top left panel shows values of 
q a ,p a , recorded when Im(z) = 0, for seven orbits; top right panel 
shows values of z when q a = for the same seven orbits (with two 
of the orbits overlapping two others, making it appear as if only 
five orbits are plotted). The fixed point in the z-plane marked 
by an 'x' is an invariant loop. When this orbit's radius is plotted 
versus its azimuth via equations (8)-(9), it gives a ID curve, as 
shown in the bottom panel. (Specifically, it gives the curve in the 
bottom panel that lies immediately above the librating regions.) 
Also shown in the bottom panel are five other invariant loops for 
five other energy values. Compare this plot with Fig. 3 obtained 
from full numerical integrations. 

orbit, as given by equations (8)-(9). For example, the 
fixed point marked with an 'x' in the upper-right panel 
of Figure 5 gives rise to one of the curves in the lower 
panel of that figure. The five other invariant loops shown 
there were found similarly, but with different values of 
the energy. 

Having described how to find the invariant loops, we 
turn now to the full 4:1 Hamiltonian (eq. [A26]). Figure 
6 shows orbits when ec = 0.01. In the sample surface 
of section shown, the fixed point is no longer at z = 
because the Lindblad resonances give a non-zero forced 
eccentricity that is in addition to the forced secular value. 
The invariant loops in the right panel of the figure are 
similar to those with only the corotation resonance (Fig. 
5), but are slightly more distorted. When Charon's ec- 
centricity is increased to ec — 0.025, the invariant loops 
become highly distorted (Figure 7). By ec = 0.03, many 
of the librating invariant loops no longer exist — the fixed 
points break up into a sea of chaos (Fig. 8). At even 
higher ec, there are no corotation librations left. 

In §3.2 we found from direct integration of Newton's 
equations that the corotation islands disappear when 
ec > 0.024. Although the critical value we find in the 
present section ~ 0.03 is similar, we speculate that the 
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Fig. 6. — Surface of section (left panel) and invariant loops 
(right panel) from integration of the full Hamiltonian (cq. [A26]), 
when ec = 0.01. The fixed point in the surface of section that is 
marked with an 'x' produces the invariant loop that is just above 
the librating retions. Notice that the inclusion of the Lindblad 
terms has shifted the value of z for the fixed point (eq. [11]). 
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Fig. 7.— Same as Fig. 6 but for e c = 0.025. The fixed point 
in the surface of section marked with an 'x' produces one of the 
librating invariant loops. New fixed points have appeared in the 
surface of section; these give invariant loops that self-intersect (not 
shown). 

reason for the discrepancy is that our numerical coeffi- 
cients are not very precise. For example, we have ne- 
glected terms that are 0{M c /M P ) ~ 0.1. 

Why does chaos appear in the Hamiltonian model of 
equation (A26) when ec > 0.03? Chaos often appears 
when separatrices of two neighbouring resonances over- 
lap (Chirikov 1979). However, it is difficult to apply this 
resonance overlap criterion to our 4:1 Hamiltonian, be- 
cause the separatrices have non-trivial shapes in four- 
dimensional phase space. Instead, we give two semi- 
quantitative explanations. 

First, we calculate the eccentricity that the first two 
Lindblad terms (ci,C2) force on a particle at the center 
of corotation. The equation of motion for z at q a = 
(center of the corotation resonance) is 

z = —2i/j, (B 2 z + cie c + c 2 ecRe(z)) , (10) 

setting c 3 = 0. The forced eccentricity is determined by 
setting z = 0, which yields 

2 

Re(z for ccd) = - „ C +° • (11) 
The forced eccentricity is infinite when ec is given by 

e c * = -— =0.031 (12) 

c 2 

We note that Ward & Canup (2006) perform a calcula- 
tion similar to our equation (11) in their Supplementary 
Notes, although they include only the first Lindblad term 
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Fig. 8. — Same as Fig. 6 but for ec = 0.03: Two left panels 
show surfaces of section at two energies. The fixed point (marked 
'x') in the upper-left panel produces the outermost librating orbit 
shown in the plot of invariant loops. When the energy is decreased 
by a small amount to 2.5 X 10~ 6 , instead of producing an invariant 
loop with larger libration, the fixed point disappears into a sea of 
chaos (lower-left panel). The chaotic orbit gives rise not only to the 
splattering of points shown near z = 0.017; its surface of section 
also has three chaotic islands at large \z\. These islands are not 
seen in the figure, because they are off the scale. 

(i.e., they effectively set c 2 — c% — 0). We also note that 
if we do not set C3 = 0, the forced eccentricity no longer 
diverges, although it still becomes large when ec ~ ec*- 
For our second explanation, we argue that chaos is 
likely to occur when the center of the second Lind- 
blad resonance (the c 2 term) overlaps the center of 
the corotation resonance in the q a ,p a plane. To find 
the center of the second Lindblad resonance, we switch 
from z to the canonical variables q e ,p e , defined via 
y/2p e e tq " = z. Then, the Hamiltonian for the second 
Lindblad resonance becomes H = — 24p^ — 2p e puB 2 — 
2/zc 2 ecP e cos(g a + 2q e ). The center of this resonance oc- 
curs where q a + 2q e = 0. Employing Hamilton's equa- 
tions, we find = q a + 2q e = — 48p a — 4/i£>2 — 4/xC2ec- 
Therefore the center of this resonance overlaps the center 
of the corotation resonance (which is at p a = 0) when ec 
is given by ec* ■ 

5. ALTERNATIVE FORMATION SCENARIOS 

How did Nix and Hydra form? Ward & Canup (2006) 
argue that these moons are byproducts of the impact that 
formed Charon. But here we argue that this is not the 
case. If Nix and Hydra were byproducts of the impact, 
one might imagine three possibilities for how they ended 
up in their current orbits, with very small eccentricities: 

• They might have been directly ejected to their cur- 
rent semimajor axes with large eccentricities, which 
were subsequently damped by tides. However, the 
timescale for them to damp their eccentricities by 
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tidal interaction with Pluto is longer than the age of 
the Solar System (Stern et al. 2006). 9 In addition, 
they likely could not have damped their eccentric- 
ity by exciting Charon's eccentricity with Charon 
in turn tidally damping its own eccentricity (Lith- 
wick & Wu 2008). 

• They might have been ejected from the Charon- 
producing impact to their current scmimajor axes 
along with many small particles. If these small 
particles formed into a collisional disk, the disk 
could have damped Nix and Hydra's eccentricities. 
However, a post-impact collisional disk probably 
could not have extended to such a large distance 
(<~ 50Rp) from Pluto: simulations of the Pluto- 
Charon impact give much more compact disks 
(Canup 2005). 

• As proposed by Ward & Canup (2006), they might 
have been damped by a collisional disk much closer 
to Pluto, and then been migrated outward by 
Charon. However, we have shown in this paper 
that both Nix and Hydra could not be resonantly 
migrated outward by Charon in the 4:1 and 6:1 
resonances, respectively. Is it possible that Charon 
was responsible for resonantly migrating only Nix 
(ec < 0.024), while Hydra was transported out- 
ward in Nix's 3:2 corotation resonance? We have 
failed to find such an orbit using our algorithm 
(§3.2) in the presence of an eccentric Charon and 
a massive Nix. But we have yet to exclude this 
possibility with more confidence. 

Taken together, the above results suggest that Nix and 
Hydra are not byproducts of the Charon- forming impact. 
In addition, Nix and Hydra could not have been cap- 
tured from the Kuiper belt, and subsequently collision- 
lessly hardened by other passing-by bodies. Although 
that mechanism might work for other Kuiper belt bina- 
ries (Goldreich et al. 2002), it would not be consistent 
with Nix and Hydra's small eccentricities and inclina- 
tions. 

Instead, we argue that Nix and Hydra formed within a 
collisional Plutocentric disk that was composed of small 
bodies captured from heliocentric orbits. It is quite plau- 
sible that at early times in the history of the Solar Sys- 
tem there were many small bodies in heliocentric orbit. 
Three out of the four largest KBOs (including Pluto) 

9 Nix could have damped its eccentricity if it is a strcngthlcss 
rubble pile and has both tidal Love number and tidal Q factor of 



have satellites that have been suggested to form out of 
impacts (Brown et al. 2006). Charon formed out of an 
impact between two Pluto-sized objects (Canup 2005). 
Such events are highly unlikely in the environment of 
the present day Kuiper belt, where the time before the 
next impact between two Pluto-sized bodies is <~ 3 x 10 13 
yr. This implies that in the past the velocity disper- 
sion within the Kuiper belt was much smaller than it is 
today, in which case gravitational focusing would have 
enhanced the collision rate between Pluto-sized bodies. 
The most likely mechanism to cool the population of Plu- 
tos is dynamical friction with a large mass of small bodies 
(Goldreich et al. 2004). Today, these bodies have suf- 
fered collisional diminuation and no direct evidence of 
their past abundance remain. However, they must have 
been present, at least at some point during the accretion 
growth period, to form the large KBOs in less than the 
age of the Solar System (Kenyon & Luu 1998). 

The small bodies collide much more frequently and can 
be collisionally accreted to the big bodies (Sari & Gol- 
dreich 2006). Collisions remove their random momenta 
but cannot remove the net angular momentum. They 
then form a disk around the big bodies. Moons, formed 
or captured by such a disk, can be migrated inward, be 
parked near resonant locations, and have small free ec- 
centricities (Shannon et al., in preparation). The moons 
are expected to be coplanar with Charon - even if its 
nascent disk started differently, such a dissipative disk 
can quickly relax to Charon's orbital plane. 

Hydra's orbit extends only to <~ 1% of Pluto's Hill 
radius. Searches by Nicholson & Gladman (2006) and 
Steffi et al. (2006) have excluded the presence of other 
comparably-sized moons in the Hill sphere. Why do 
present moons occupy such a small fraction of the avail- 
able space? We suspect the moon-harboring disk may 
be limited in size, determined by the net angular mo- 
mentum of the accreted small bodies. Depending on the 
velocity dispersion of the small bodies in the circumso- 
lar disk, the size of the accreted circum-Pluto disk might 
be significantly smaller than Pluto's Hill radius. Further 
exploration is underway. 

Pluto is not special. Similar accretion disks may also 
have arisen around other large Kuiper belt objects. Do 
they also possess small moons? Hopefully the moons of 
Kuiper belt objects can teach us about the early history 
of our Solar System, when the planets — and the KBO's 
themselves — were being formed. 

order unity; but Hydra could not have, even in this extreme case. 



APPENDIX 

A. HAMILTONIAN NEAR THE 4:1 RESONANCE 

Consider a masslcss test particle that is coplanar with the Pluto-Charon binary, near its exterior 4:1 resonance. Pluto 
and Charon's mutual orbit has orbital elements {ac 1 ec, Ac, ^c, n c}, with Q(Mp+Mc) — n c a^ and Ac = const+npt. 
We set txjc — without loss of generality. Since Mc <C Mp, the effect of Charon on the particle may be treated as 
a perturbation to the effect of Pluto. The particle has Pluto-centric orbital parameters 10 {a, e, A, zu}, and its energy 
per unit mass is 

E= _GMp_GMc n 

2a a x ' 

Throughout most of this paper, we use Jacobi co-ordinates, where the test particle's orbital elements arc relative to the barycenter of 
Pluto and Charon. But in this Appendix and in §4, we employ Pluto-centric elements because this is traditional when using a disturbing 
function (Murray & Dcrmott 1999)-even though it would be simple to use Jacobi elements instead (Lithwick & Wu 2008). 
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TABLE 2 

Leading cosine terms in TZ (eq. [Al]) near. 4:1 resonance, from Appendix B in Murray & Dermott (1999) 



Cosine Argument 


Coefficient 




ra 


Bie c e 


Bi 


= -0.0902 





B 2 e 2 


B 2 


= 0.0898 


4A- A c 


C e3 


Co 


= -0.285 


4A - Ac - to 


Cie^e 


Ci 


= 1.82 


4A - Ac - 2ro 


C 2 ece 2 


c 2 


= -3.82 


4A - A c - 3to 


C 3 e 3 


c 3 


= 0.640 



where TZ can be expanded as a sum of cosine terms. In Table 2 we list the coefficients and arguments of the leading cosine 
terms near the 4:1 resonance. To translate our notation to that of Appendix B in Murray & Dermott (1999), from which 
we extracted the numerical constants in the table, TZ = TZd + a~ 2 TZi, where a = ac/a, and {Bi, B2, Co, C\, C2, C3} = 
{/io, /2j /s2> /s3> /s4> /ss — VC^c* 2 )}. For the numerical constants, we set a = ac/a rcs , where ct ros is the semimajor 
axis at nominal 4:1 resonance (eq. [A9]). 

The Hamiltonian is equal to the energy, after replacing the particle's orbital elements with canonical variables 
H(A, A, r, 7) = E. We adopt Poincare canonical variables {A, A, T, 7}, where 

A=(gM P a) 1/2 (A2) 

T=(gM P a) 1/2 e 2 /2 (A3) 
7=-w, (A4) 

(Murray & Dermott 1999), dropping terms 0(e 4 ) from T. 

To simplify the Hamiltonian, we employ a number of variable transformations (Holman & Murray 1996; Murray & 
Dermott 1999). First, we absorb the time-dependent parameter Ac = npt+const into A and shift A so that it vanishes 
at the nominal 4:1 resonance, 

{A',A'}^|^^,4A-A C | , (A5) 

where A rcs is a constant to be determined. Since the generating function for this transformation is F — (A' + 
A res /4)(4A — Ac), the new Hamiltonian is H + d t F, 

H(A ' A ' r ' 7) - "2(4A' + A rcs )2 " (A + T" )HC (A6) 

setting a to its value at nominal resonance, a — a ICS (eq. [A9]), everywhere except in the first two terms. To choose A ros , 
we require that A' = at the nominal 4:1 resonance, which occurs where = {d/dt){A\ — \ c ) = dX'/dt = dH/dh! . 
Since 

dH _ 4G 2 Mp 

dK' A '=o ~ a? cs 

we set 



n c -» , (A7) 



A rcs = (^Mp/nc) 1 / 3 . (A8) 
The value of a at nominal resonance is from equation (A2), 

°'- = 42/ ^(iOT) V> ' (A9) 

Expanding the Hamiltonian to second order in A' and dropping the constant term 

&res ^*- r es ^res 

We may rescale the momenta and Hamiltonian by the same constant factor without altering the equations of motion. 
Rescaling by A rcs , the Hamiltonian becomes 

2 

H(p a ,q a ; e ^,-zv) = - r f(24pl + »1l) (All) 

l"^ (A13) 

9a = A' = 4A-A c (A14) 



where 
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The canonical momenta are p a and e 2 /2, and their corresponding conjugate co-ordinates are q a and —vj. 

Next, we transform {e 2 /2, — vo} to remove the forced secular eccentricity. It simplifies the algebra to switch to the 
complex canonical variable (Strocchi 1966; Ogilvie 2007) 

Z = ee~ m , (A15) 

which is the usual complex eccentricity. Hamilton's equations for {e 2 /2, — w} are now expressed as 

dZ n dH 

Tt= 2 ^ (A16) 

where Z* is the complex conjugate of Z. The secular part of 1Z in Table 2 is 

TZ scc ^e c B 1 ^-^ + B 2 \Z\ 2 (A17) 



Bo 



B 1 2 



YeC 2B 2 

after dropping a constant. Therefore, we transform to the variable 



(A18) 



Z + e c~, (A19) 



2 

which is the (complex) free eccentricity; the constant offset is the forced eccentricity, ef orce d = —ec{B\/2B 2 ). Hamil- 
ton's equation for z is clearly the same as for Z (eq. [A16]), i.e., the transformation is canonical. 
Under transformation (A19), the resonant part of 1Z becomes 

ll ICS =Rc (e^ (C e 3 c + de 2 c Z + C 2 e c Z 2 + C 3 Z 3 )) (A20) 

= Rc(e tqa (c e 3 c + c^z + c 2 e c z 2 + c 3 z 3 )) (A21) 

where, defining [5 = —B\/2B 2 , 

c = C + Crf + C 2 f3 2 + C 3 (3 3 = -0.26 (A22) 

c 1 = Ci + 2C 2 /3 + 3C 3 /3 2 = -1.5 (A23) 

c 2 = C 2 + 3C 3 f3 = -2.9 (A24) 

c 3 = C 3 = 0.64 (A25) 

Collecting results, the Hamiltonian is 

H( Pa ,q a ;z) = -24p 2 -/iB 2 |z| 2 -^Re(e l9 » (c a e 3 c + Cl e 2 c z + c 2 e c z 2 + c 3 z 3 )) , (A26) 

after dropping the prefactor nc/4, which means that time is now measured in units of 4/nc- Hamilton's equations of 
motion are 

dH dH . n dH fl . 

P * = Wa' qa = -W» ; (A27) 
B. TIDAL DISSIPATION 

Assuming orbits and spins are coplanar, tidal dissipation due to tides raised on Charon by Pluto are described by 
the following evolution equations (Hut 1981) 



6&2C s { Rc\ a ( r /-, 2\3/2j-^ 



+ «) ^ 77-^™ h - (1 - e 2 )3/ 2 / 2 ^ , (Bl) 



c Tc \ a J (1 - e 2)15/2 

c = -^ (1 + nvJ (l-e 2 )W 2 p-18 (1 ~ e) /4 V^' ,L ' 2) 



da 
~dt 

de 
dt _ 

dO c _3fc 2C g 2 n ^c\ 6 f . n 2x3/2 , fi c\ 

^-^^(T^f^J V 2 - (1 - e) /5 V/' (B3) 

where q = Mp/Mc is the mass ratio, f^c is Charon's spin rate, n = yjG{Mp + Mc)/a 3 is the orbital mean motion, 
k 2 c is the tidal Love number of Charon, and r g c its radius of gyration related to the moment of inertia by / = r 2 mR 2 , 
with r g ~ 0.6 for a uniform density sphere. The values for /j are all of order unity, /j = 1 + 0(e 2 ) + C(e 4 )..., 
with their exact expressions listed in Hut (1981). In this tidal model, Tc — R c /GMcTc where tq is the (assumed) 
constant tidal lag time. To connect this model to that of Goldreich & Soter (1966) which assumes a constant lag phase 
(and a constant tidal Q factor), we let the tidal lag time tc = l/(2nQc), with Qc being Charon's tidal dissipation 
quality factor. This choice is somewhat arbitrary: one may also assume, e.g., tq — l/[2(f2c — n )Qc\- The tide raised 
by Charon on Pluto causes similar effects but with all subscripts C substituted by P and the mass ratio inverted, 
q = Mc/Mp. The net orbital evolution is given by the sum of both tides. 
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The Love number for a uniform density sphere is (see, e.g., Dobrovolskis et al. 1997; Murray & Dermott 1999) 

*■ - r& < B4 > 

where jl = l9/i/(2pgR) is the effective rigidity of the body, p the material strength, p the density and g the gravitational 
acceleration. If the body is a rubble pile, jl <C 1 and k 2 ~ 3/2. However, if the body is a consolidated ice sphere with 
/i~4x 10 10 dync/cm 2 ,'/c 2 c ~ 0.005 and k 2 p ~ 0.05 (Dobrovolskis et al. 1997). 
The following are numerical estimates for the various timescales, t x = x/\x\, 



'alp ' 



Qc\ ( k 2 c \ / 



6.5 



r.| o ~3xl0V.^)^j (j^j , (B7) 

^*^<M)(^y\^T- (B8 > 

^-•"M^X^rGsb)" (B9) 



/ \ 6.5 

f ac 



T „ p . 6xl0Vs ^J^J ^_±_J , (B10) 

where we have scaled Qp and Qc by values typical of solid bodies and have taken the limit e«l. 

Charon's spin is quickly (pseudo-)synchronized with the orbit, and Pluto's spin synchronization is on-going until the 
system reaches the final state of double synchronization. So most of the a evolution is contributed by tides on Pluto 
(eq. B6). 

During the orbital expansion, the tide on Pluto increases e while the tide on Charon decreases it. With our choice 
of parameters, the former dominates over the latter. Once Charon reaches its equilibrium location, both tides damp 
the eccentricity. 

If the spin direction of body is misaligned with the orbit normal, it will be tilted to alignment in roughly the spin 
synchronization time (Hut 1981). 
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